# -*- coding: utf-8 -*-
#优化统计小区人口

import utm
import os
import pyproj
from itertools import islice
import numpy as np
from osgeo import gdal
import pandas as pd
from shapely.geometry import Point
from geopandas import GeoDataFrame
from shapely.ops import transform
from shapely.geometry import Polygon

toponode={}
lat,lon = [],[]
dataSet = []

#栅格处理函数，输入矩形经纬度坐标
#返回zone_point（经纬度坐标）
def raster_operation(xmin, ymax, xmax, ymin):
    # 读取部分数据集
    # print('******************')
    # print('开始读取区域栅格...')
    # os.chdir(r'E:\MyWork\0_SD_Local\project\03_cq\data')
    in_ds = gdal.Open(pop_lig_file)
    in_gt = in_ds.GetGeoTransform()  # 获得输入栅格数据的地理变换（geotransform）

    inv_gt = gdal.InvGeoTransform(in_gt)  # 获得地理变换（geotransform)的逆
    if gdal.VersionInfo()[0] == '1':
        if inv_gt[0] == 1:
            inv_gt = inv_gt[1]
        else:
            raise RuntimeError('Inverse georansform failed')
    elif inv_gt is None:
        raise RuntimeError('Inverse geotransform failed')

    offsets_ul = gdal.ApplyGeoTransform(inv_gt, xmin, ymax)  # 获得左上角和右下角偏移坐标（图像坐标）
    offsets_lr = gdal.ApplyGeoTransform(inv_gt, xmax, ymin)
    off_ulx, off_uly = map(int, offsets_ul)
    off_lrx, off_lry = map(int, offsets_lr)
    # print('左上角图像坐标:')
    # print(off_ulx, off_uly)
    # print('右下角图像坐标:')
    # print(off_lrx, off_lry)

    rows = off_lry - (off_uly - 1)  # 获得要提取区域的行数和列数
    colums = off_lrx - (off_ulx - 1)  # 图上坐标中左上角的x,y坐标最小
    # print('行数：' + str(rows), '列数：' + str(colums))

    in_band = in_ds.GetRasterBand(1)
    data = in_band.ReadAsArray(off_ulx, off_uly, colums, rows)  # 研究区域的栅格数据，Numpy二维数组格式
    # print("区域栅格数组：\n"+str(data))

    # 研究区域的栅格数据处理
    # print('开始栅格数据处理...')
    xyCoors = []
    GRID_CODE = []
    # with rasterio.open('E:\MyWork\\0_SD_Local\project\\03_cq\data\china.tif') as src:
    #     # print(src.xy(off_uly,off_ulx),src.xy(off_lry,off_lrx))
    #     for r in range(off_uly, off_lry + 1):  # 获得区域栅格的中心点坐标以对应的值
    #         for c in range(off_ulx, off_lrx + 1):
    #             xyCoors.append(src.xy(r, c))
    #             GRID_CODE.append(data[r - off_uly, c - off_ulx])

    for r in range(off_uly, off_lry + 1):  # 获得区域栅格的中心点坐标以对应的值
        for c in range(off_ulx, off_lrx + 1):
            # xyCoors.append(src.xy(r, c))
            xyCoors.append(gdal.ApplyGeoTransform(in_gt, c, r))
            GRID_CODE.append(data[r - off_uly, c - off_ulx])


    # print('坐标点个数：')
    # print(len(xyCoors))
    # print(xyCoors)
    # print('点的值个数：')
    # print(len(GRID_CODE))
    # print(GRID_CODE)

    s1 = pd.Series(xyCoors)  # 创建DataFram
    s1.name = 'geometry'
    s2 = pd.Series(GRID_CODE)
    s2.name = 'GRID_CODE'
    df = pd.DataFrame([s1, s2], index=['geometry', 'GRID_CODE'])
    df = df.T

    geometry = [Point(xy) for xy in df['geometry']]  # DataFram转化为GeoDataFram
    df = df.drop(['geometry'], axis=1)  # 从行或列中删除指定的标签。
    crs = {'init': 'epsg:4326'}
    zone_point = GeoDataFrame(df, crs=crs, geometry=geometry)

    return zone_point

#经纬度转UTM，输入（经度，纬度）
def wgs2proj(y_84,x_84):
    (x_coor, y_coor, temp_1, temp_2) = utm.from_latlon(x_84, y_84)
    if not (zone_number == temp_1):                                                 # 存在跨zone的点
        p1 = pyproj.Proj(proj='utm', zone=temp_1, datum='WGS84')
        p2 = pyproj.Proj(proj="utm", zone=zone_number, datum='WGS84')
        x_coor, y_coor = pyproj.transform(p1, p2, x_coor, y_coor)                  #转换带号
    return x_coor, y_coor

#文件夹目录处理
f = open("../files/runinfo_vc6.txt")
next(f)
for d in f:
    d = d
zoneDic = str(d)+"/data/zone.txt"
toponodeDic = str(d)+"/data/toponode.txt"
latlon_utmDic = str(d)+"/data/latlon-utm.txt"
TripPopulationDic = str(d)+"/data/TripPopulation.txt"
autoZoomParaDic = str(d)+"/cache/autoZoomPara.txt"

#灯光人口数据文件路径
f = open(autoZoomParaDic)
next(f)
d2 = f.readline()
d3 = f.readline()
pop_lig_type = d2.split()[1]
pop_lig_file = d3.split()[1]

f = open(latlon_utmDic)                   #获得zone和band
d1 = f.readline()
zone_number = int(d1.split()[5])
zone_letter = d1.split()[6]

for line in open(toponodeDic):
    toponode[line.split()[0]]=(float(line.split()[1]),float(line.split()[2]))     #将拓扑点的坐标做成一个字典，键为ID，值为坐标的元组

print("--开始遍历小区--")
for line in islice(open(zoneDic),1,None):
    temp=line.split()
    p=map(lambda x:toponode[x],temp[1:])                                                   #用MAP函数对所有拓扑点ID找到其对应的UTM坐标
    tz = Polygon(list(p))                                                                  #将连续的坐标点转换成面要素

    q = map(lambda x: toponode[x], temp[1:])
    myarray = np.asarray(list(q))
    coor_y = myarray[:,0]
    coor_x = myarray[:,1]
    codeNum = len(myarray)                                                                  #单个小区边界中点的个数

    # UTM转经纬度
    for i in range(0, codeNum):
        lat_temp, lon_temp = utm.to_latlon(coor_y[i], coor_x[i], zone_number, zone_letter)
        lat.append(lat_temp)
        lon.append(lon_temp)
    xmin, ymax = min(lon), max(lat)                                                         # 矩形左上角经纬度坐标
    xmax, ymin = max(lon), min(lat)                                                         # 矩形右下角经纬度坐标

    # 调用栅格数据处理函数
    zone_point = raster_operation(xmin, ymax, xmax, ymin)
    # 经纬度转UTM
    zone_point.geometry = zone_point.centroid.apply(lambda z: transform(lambda x, y: wgs2proj(x, y), z))

    # print(tz)
    # print(zone_point)
    inzone=zone_point[zone_point.within(tz)]                                                #判断在小区内部的点
    # print(inzone)
    a=sum(inzone.GRID_CODE)
    dataSet.append([int(temp[0]),a])
    print(str(temp[0])+"小区人口："+str(a))


print("--小区数据统计完成--")
with open(TripPopulationDic, 'w') as f:
    f.write('小区编号 小区常住人口\n')
    for i in range(len(dataSet)):
        for j in range(len(dataSet[i])):
            f.write(str(dataSet[i][j]))
            f.write(' ')
        f.write('\n')
f.close()
print("--TriPopulation.txt已生成--")



